Hybrid method for simulating front propagation in reaction-diffusion systems 
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We study the propagation of pulled fronts in the A «-> A + A microscopic reaction-diffusion process 
using Monte Carlo (MC) simulations. In the mean field approximation the process is described 
by the deterministic Fisher-Kolmogorov-Petrovsky-Piscounov (FKPP) equation. In particular we 
concentrate on the corrections to the deterministic behavior due to the number of particles per site 
fl. By means of a new hybrid simulation scheme, we manage to reach large macroscopic values of Q 
which allows us to show the importance in the dynamics of microscopic pulled fronts of the interplay 
of microscopic fluctuations and their macroscopic relaxation. 
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When describing systems at much larger scales than 
correlation length, internal fluctuations due to the intrin- 
sic discreteness of the particles can be neglected, since 
they only account for a correction typically of the order 
of ft -1 / 2 , where ft is the number of particles in a corre- 
lated volume. However, in some situations the dynamics 
of the system spans different scales which give rise to 
a strong dependence of its macroscopic features on the 
microscopic details of its constituents, even in the limit 
ft — > oo. Relevant instances of this phenomena are the 
dynamic contact angle problem Q , evolution of a fracture 
tip , dendritic growth Q , and the flow of a gas through 
a microscopic channel Q . In this paper we highlight and 
study another important example namely, the effect of in- 
ternal fluctuations in the macroscopic dynamics of pulled 
fronts pj, . Specifically we consider the propagation of 
pulled fronts in reaction-diffusion microscopic problems, 
like the A <-> A + A scheme 0,110 Continuum de- 
scription of the system is possible in the reaction-limited 
regime where ft is large enough so the reaction is well 
stirred within each correlated volume • I n this case, 
the density p{x, t) of particles per correlated volume is 
described, in the limit ft — > oo, by the FKPP equation 

M 



dp n d 2 p 



Ot 



k 2 p 2 . 



(1) 



This equation has travelling-wave solutions of the form 
p = p(x — vt) which invade the unstable phase p[pd) = 
from the stable phase p{— oo) = ki/k^ and travel with 
velocity v > vq = 2V Dk\. Of particular interest is the 
solution with velocity Vq, since it is dynamically selected 
for a broad class of initial conditions. Moreover, i> is 
the linear spreading speed of infinitesimal perturbations 
around the unstable state. Thus, fronts with velocity vq 
are essentially "pulled along" by the growth and spread- 
ing of small perturbations in the leading edge x 3> vt 
where p < 1. This sensitivity causes also the absence 



of a typical macroscopic length and time scale in which 
perturbations around the asymptotic solution with ve- 
locity vq are damped For example, the velocity of 
the front starting from an steep enough initial condition 
approaches the asymptotic value like a power law 



v(t) = V 



3 
2q^t 



+ o(t- 3 / 2 ) 



(2) 



where qo = vq/2D. 

In particle models however, the continuum description 
given by the FKPP equation breaks down at p ~ I/O 
where internal fluctuations are important. Since pulled 
front dynamics are sensitive to infinitesimal events at 
p <C 1 we expect macroscopic properties to depend 
strongly on f2 when il — > oo. For example, by neglecting 
microscopic fluctuations and mimicking the discreteness 
of particles by imposing an effective cutoff in the FKPP 
equation at p = O -1 , Brunet and Derrida [ill Il2l| ob- 
tained that the velocity is given by 



v(oo) = Vq = v ~ 



vqICv 

inn 
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where JC V is a constant. As expected the correction to 
the macroscopic velocity of the front is very strong: in a 
macroscopic volume of 10 23 particles, it is still of 0.3%. 
Combining Eqs. J2J and © we can easily infer that the 
typical time scale of microscopic pulled fronts is given by 
the condition v(tq) ~ vq, that is 



In 2 ft 



(4) 
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Thus, pulled fronts in reaction diffusion particle models 
do have a typical time scale as opposite to those of the 
FKPP equation, although it is set by microscopic details 
and diverges in the limit ft — > oo |14| . 

Numerical confirmation of these predictions in general 
reaction-diffusion particle models is difficult, since the 
observation of the functional dependence in @ requires 
typical simulations up to In ft ~ 10 2 which are compu- 
tationally unfeasible. However, for a particular model in 
which particles undergo non-local diffusion movements, 
Brunet and Derrida where able to perform simulations 
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up to ft ~ 10 150 and check with high accuracy the pre- 
diction (0) 0,^3- Moreover, they found also that the 
front diffuses in time and that the diffusion coefficient 
behaves as ^1 



ICi 



In' 5 ft 
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where JCd is a constant. While the velocity correction © 
can be easily understood in terms of an effective cutoff 
in the FKPP equation , and simulations for moderate 
numbers of particles (ft < 10 10 ) in reaction-diffusion par- 
ticle models seem to be compatible with © 6.l^.ll3l lT5| . 
the functional dependence of Z?q has only been observed 
in the non-local model of Brunet and Derrida. Although 
there is an heuristic argument for a specific model to get 
the In -3 ft dependence 0, the situation clearly remains 
unsatisfactory, since there is only one empirical observa- 
tion of J5J|. Thus, our purpose in this paper is to simulate 
the A «-> A + A model for very large number of particles 
in order to check both © and JSJ and get some insight 
into the dynamics of pulled fronts in microscopic particle 
models. 

The A <-> A + A model in one dimension consists of 
particles on a lattice with spacing Ax in which the num- 
ber of particles at site i, Ni(t) is unbounded |g. Reaction 
events take place on-site, while diffusion drives particles 
to nearest neighbors positions. Particles annihilate with 
rate tr, give birth to another one with rate 7 and dif- 
fuse with rate D. In equilibrium, the average number 
of particles per site is ft = 7/cr, and when ft — ► 00 the 
system is described by Eq. QJ with p(x,t) — N(x,t)/Q 
and ki = k-2 = 7. Early MC simulations of this model 
[j| showed that indeed when ft 3> 1, FKPP pulled fronts 
emerge. Specifically, it was found that both the velocity 
correction and the diffusion coefficient of the front decay 
like ft -1 / 3 for ft < 10 6 , a scaling which has been observed 
in other models for moderate values of ft [lq. However, 
in order to observe the scalings © and (JSJ much larger 
numbers of particles are needed (typically ft ^> fO 10 ) 
which can not be attained in standard MC simulations. 

To reach larger number of particles in the A <-> A + A 
model we note that, since the dynamics of pulled fronts 
are very sensitive to the dynamics of the system close 
to the unstable state p ~ 0, a correct description of 
the reaction-diffusion microscopic problem is only needed 
there, where fluctuations are important |T^. Away from 
the unstable state, i.e. when p = 0(1), the number of 
particles is big enough so that fluctuations are negligible 
and the system can be safely described by macroscopic 
descriptions like JTJ- Thus we propose to split the dy- 
namics of the microscopic model into two different de- 
scriptions: given a mesoscopic number of particles N* 
with ft ^> N* ^S> 1, at any time step t we identify the po- 
sition i* as the smallest value of i for which ATj(i) < N* . 
In the region in which fluctuations can be safely ignored, 
that is when i < i* we update the number of particles 
using a numerical approximation of Eq. , while we use 
MC methods in the region i > i*. To complete the al- 




FIG. 1: Snapshot of the front profile for one realization of 
the hybrid method with Q, = 10 20 (upper panel). Below: 
schematic view of the front Ni(t) close to i*. Open symbols 
are in the macroscopic region and full symbols in the micro- 
scopic region. For each point lattice i we define the flux of 
outgoing F^(t) and incoming F~ (t) particles. In the figure 
only those for i* are shown. Lines are guides to the eye. 



gorithm, boundary conditions at i = i* should be given. 
Since only diffusion couples the dynamics between differ- 
ent sites, we implement the boundary condition through 
the conservation of fluxes of particles through the bound- 
ary, similarly to other MC hybrid methods 0, 0, 0] . 

To this end, if we define at each site i the flux of in- 
coming particles F~ (t) and of outgoing particles (t) 
(see Fig. , the Euler approximation with time step At 
of Eq. QJ reads 



Ni{t + At)-Ni(t) _ Fr(t)-F+(t) 



At 



Ax 



with 



+ lNi (t)-aN?(t). 

(6) 



^[Ni-i(t) - Ni(t)] 
£-[Ni(t) - N i+1 (t)} 



(7) 



Obviously, conservation of the number of particles re- 
quires that Fr(t) = F+^t) and F+(t) = F~ +l {t). Thus, 
our algorithm evolves as follows: for a given mesoscopic 
time step At, we update the microscopic region (i > i*) 
using time continuous MC until the time of the simu- 
lation At exceeds At. Note that the typical time step in 
the MC simulation is given by St^ 1 ~ A*[l + log(ft/A*)] 
which is smaller than At in our simulations. Thus, 
several MC events take place until the MC simulation 
time At exceeds At, which makes the real At different 
for any time step. Any MC event in which a particle 
jumps into the macroscopic region (i < £*) is recorded 
in the variable N~ and the particle is removed from 
the MC simulation. Wc then update sites i < i* — 2 
using Eq. © and Q. Finally, the number of parti- 
cles at the boundary site — 1 is updated using equa- 
tion JfJJl but with F 1 i f_ 1 calculated according to the MC 
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FIG. 2: Asymptotic velocity correction as a function of the 
logarithm of the number of particles Q. Open symbols are 
results for the full MC simulation while full symbols are for 
the hybrid scheme. Dashed line corresponds to the scaling 
given by with K. v given by (|5J. Error bars are not shown 
when smaller than the symbol size. Inset: time evolution 
of the instantaneous velocity of the front (solid lines) with 
n = 10 20 , 10 30 and 10 40 from top to bottom. Dashed line is 
the prediction given by Eq. @. 



recorded number of jumps, N . Specifically we take 
F+_ x = (D/Ax)Ni*-x{t) - N~/(AtAx). Since we 
should get that F^{t) = F^_ x {t), we update the number 
of particles at site i* to fulfill this condition on average: 
Ni,{t + M) = N i .(t + M)+TL A tD Nt ._ 1 tt)/* ai where II A is 
a Poisson random number with mean A. This completes 
a time step At in the algorithm. 

The condition for a sharp interface between the macro- 
scopic and the microscopic region at i* can be relaxed by 
introducing a buffer region 16, 17] . Moreover, fluctu- 
ations can be also considered in the macroscopic region 
by adding an internal noise source to the FKPP equation 
[pa.l20| . However, for large enough TV* and small enough 
At our results do not differ from those of these algo- 
rithm refinements. In our simulations we take At = 10~ 4 , 
N* = min{10 4 , Cl/2} and D = 7 = Ax = 1. 

Results for this hybrid scheme are shown in figures 
121 and |21 and compared with full MC simulations up to 
CI = 10 5 . Both the correction to the velocity of the front 
© and the diffusion coefficient |J5J agree, up to statisti- 
cal fluctuations, with those of the full MC simulations, 
which supports the validity of our algorithm. Note that 
when CI — > 00 our algorithm reduces to the Euler approx- 
imation of the FKPP equation. Thus, the velocity vq in 
© is given by the solutions of the equations [2lJ 

Voe - qo v At = „2sinh 9o , 

e - qo v Ai = 2At[(cosh 9o -l)-l], (8) 

and the velocity correction coefficient obtained using the 
efficient deterministic cutoff argument of is given by 

JC V = ir 2 q (e VoqoAt cosh q - V%At/2). (9) 




log 10 (flD 

FIG. 3: Diffusion coefficient of the front Df as a function of 
the logarithm of the number of particles Q. Symbols are as in 
figure Dashed line is a fit of the last points to the scaling 
form © with K,d = 26.5. Dotted-dashed line is the scaling 
Df ~ cr -0 ' 32 of 8j. Error bars are not shown when smaller 
than the symbol size. Inset: Velocity as a function of time for 
one realization of the algorithm with Q = 10 20 . 



In the limit CI — > 00 we observe in Fig. [21 that our re- 
sults tend to the scaling JSJ together with the solutions 
© and 10. However, a strong deviation of our results 
for the predicted scaling J2J) even for large values of 
is observed, a fact that also was present in the Brunct 
and Derrida model |12j . Asymptotic convergence of mi- 
croscopic pulled fronts in the A «-> A + A towards the 
solution of the FKPP equation is also observed in the 
inset of Fig. [21 in which we plot the time dependence of 
the velocity for different values of CI. As expected, in our 
simulations the velocity decays accurately like (JSJ until 
it saturates to a constant value given by 1(5} . 

Regarding the diffusion coefficient, our results for the 
A <-» A + A confirm the scaling ||SJ found in 12] . Note 
that the results for small values of £1 agree with the 
scaling ~ il -1 / 3 found in the initial studies of the 
A «-> A + A model [8]. To understand the origin of 
the functional dependence of the diffusion coefficient, we 
show in the inset of figure [3| a typical realization of the 
instantaneous velocity of the front as a function of time. 
As we can see, long-lived fluctuations take place at the 
front, whose origin is in the large relaxation time of mi- 
croscopic pulled front dynamics, tq ~ In 2 CI. In order to 
check this, we have measured the time correlation of the 
instantaneous velocity of the front for different values of 
CI: 

C v (t) = ([v(s + 1) - v n ][v(s) - v n \). (10) 

Our data (see figure 0} indicate that the velocity corre- 
lation scales like 
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FIG. 4: Velocity time correlation as a function of time for 
different values of the number of particles Q. Inset shows the 
scaling given by Eq. JTT1 . 



where G(x) is a scaling function. As expected, the relax- 
ation time of velocity fluctuations is given by to and, in 
particular, Eq. is consistent with the scaling of Dq 
given by © using the Kubo formula 

Do ~ lim / + C v (t')dt' ~ (12) 

The observed scaling (jl 1|) gives us some insight about 
the effect of internal fluctuations in microscopic pulled 
fronts: at small densities, fluctuations in the number of 
particles N{(t) become important (for example, note in 
Fig. the presence of particles well ahead of the tip of 
the front). Equation suggests that those fluctua- 
tions have a strength proportional to 1/ In 5 Q. In the ex- 
istence of a typical macroscopic scale, those fluctuations 



would be damped almost instantaneously and the diffu- 
sion coefficient would have been proportional to In 5 fl. 
In fact, a similar result is obtained analytically for the 
coarse-grained continuous model (the stochastic FKPP 
equation |2(j) of the A <-> A + A model when standard 
perturbation techniques (which rely in the existence of 
a macroscopic time scale) are used However, micro- 
scopic pulled fronts do not have this macroscopic time 
scale and fluctuations are accommodated by the dynam- 
ics in a much larger time scale, to- The interplay between 
the microscopic fluctuations of strength In -5 and the 
time scale of order In 2 Q in which they relax is what pro- 
duces the dependence with In f2 observed in the diffusion 
coefficient. 

In summary, we have presented a new hybrid method 
for studying the dynamics of fronts in particle reaction- 
diffusion problems. This hybrid scheme allow us to 
investigate the asymptotic convergence of those micro- 
scopic models to the macroscopic description given by 
the FKPP equation (JIJ. In particular, we reproduced the 
scaling of the velocity correction with the number of par- 
ticles given by observed in [Tl| and inferred in other 
works. More interestingly, we have confirmed the pro- 
posed scaling for the diffusion coefficient JSJ) and showed 
that its origin is in the interplay of the typical relaxation 
time of microscopic pulled fronts and the strength of the 
microscopic fluctuations at small densities. 
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